{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "view-in-github"
   },
   "source": [
    "<a href=\"https://colab.research.google.com/github/NeuromatchAcademy/course-content/blob/master/tutorials/W3D1_RealNeurons/W3D1_Tutorial1.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "# Neuromatch Academy: Week 2, Day 3, Tutorial 1\n",
    "# Real Neurons: The Leaky Integrate-and-Fire (LIF) Neuron Model\n",
    "__Content creators:__ Qinglong Gu, Songtin Li, John Murray, Richard Naud, Arvind Kumar\n",
    "\n",
    "__Content reviewers:__  Maryam Vaziri-Pashkam, Ella Batty, Lorenzo Fontolan, Richard Gao, Matthew Krause, Spiros Chavlis, Michael Waskom"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Tutorial Objectives\n",
    "This is Tutorial 1 of a series on implementing realistic neuron models. In this tutorial, we will build up a leaky integrate-and-fire (LIF) neuron model and study its dynamics in response to various types of inputs. In particular, we are going to write a few lines of code to:  \n",
    "\n",
    "- simulate the LIF neuron model\n",
    "\n",
    "- drive the LIF neuron with external inputs, such as direct currents, Gaussian white noise, and Poisson spike trains, etc.\n",
    "\n",
    "- study how different inputs affect the  LIF neuron's output (firing rate and spike time irregularity)\n",
    "\n",
    "Here, we will especially emphasize identifying conditions (input statistics) under which a neuron can spike at low firing rates and in an irregular manner. The reason for focusing on this is that in most cases, neocortical neurons spike in an irregular manner."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "both",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:41.879305Z",
     "iopub.status.busy": "2021-05-25T01:12:41.878801Z",
     "iopub.status.idle": "2021-05-25T01:12:42.189089Z",
     "shell.execute_reply": "2021-05-25T01:12:42.187958Z"
    }
   },
   "outputs": [],
   "source": [
    "# Imports\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:42.195669Z",
     "iopub.status.busy": "2021-05-25T01:12:42.193773Z",
     "iopub.status.idle": "2021-05-25T01:12:42.281413Z",
     "shell.execute_reply": "2021-05-25T01:12:42.280597Z"
    }
   },
   "outputs": [],
   "source": [
    "# @title Figure Settings\n",
    "import ipywidgets as widgets  # interactive display\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "# use NMA plot style\n",
    "plt.style.use(\"https://raw.githubusercontent.com/NeuromatchAcademy/course-content/master/nma.mplstyle\")\n",
    "my_layout = widgets.Layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:42.295424Z",
     "iopub.status.busy": "2021-05-25T01:12:42.289192Z",
     "iopub.status.idle": "2021-05-25T01:12:42.298427Z",
     "shell.execute_reply": "2021-05-25T01:12:42.297901Z"
    }
   },
   "outputs": [],
   "source": [
    "# @title Helper functions\n",
    "\n",
    "\n",
    "def plot_volt_trace(pars, v, sp):\n",
    "  \"\"\"\n",
    "  Plot trajetory of membrane potential for a single neuron\n",
    "\n",
    "  Expects:\n",
    "  pars   : parameter dictionary\n",
    "  v      : volt trajetory\n",
    "  sp     : spike train\n",
    "\n",
    "  Returns:\n",
    "  figure of the membrane potential trajetory for a single neuron\n",
    "  \"\"\"\n",
    "\n",
    "  V_th = pars['V_th']\n",
    "  dt, range_t = pars['dt'], pars['range_t']\n",
    "  if sp.size:\n",
    "    sp_num = (sp / dt).astype(int) - 1\n",
    "    v[sp_num] += 20  # draw nicer spikes\n",
    "\n",
    "  plt.plot(pars['range_t'], v, 'b')\n",
    "  plt.axhline(V_th, 0, 1, color='k', ls='--')\n",
    "  plt.xlabel('Time (ms)')\n",
    "  plt.ylabel('V (mV)')\n",
    "  plt.legend(['Membrane\\npotential', r'Threshold V$_{\\mathrm{th}}$'],\n",
    "             loc=[1.05, 0.75])\n",
    "  plt.ylim([-80, -40])\n",
    "\n",
    "\n",
    "def plot_GWN(pars, I_GWN):\n",
    "  \"\"\"\n",
    "  Args:\n",
    "    pars  : parameter dictionary\n",
    "    I_GWN : Gaussian white noise input\n",
    "\n",
    "  Returns:\n",
    "    figure of the gaussian white noise input\n",
    "  \"\"\"\n",
    "\n",
    "  plt.figure(figsize=(12, 4))\n",
    "  plt.subplot(121)\n",
    "  plt.plot(pars['range_t'][::3], I_GWN[::3], 'b')\n",
    "  plt.xlabel('Time (ms)')\n",
    "  plt.ylabel(r'$I_{GWN}$ (pA)')\n",
    "  plt.subplot(122)\n",
    "  plot_volt_trace(pars, v, sp)\n",
    "  plt.tight_layout()\n",
    "\n",
    "\n",
    "def my_hists(isi1, isi2, cv1, cv2, sigma1, sigma2):\n",
    "  \"\"\"\n",
    "  Args:\n",
    "    isi1 : vector with inter-spike intervals\n",
    "    isi2 : vector with inter-spike intervals\n",
    "    cv1  : coefficient of variation for isi1\n",
    "    cv2  : coefficient of variation for isi2\n",
    "\n",
    "  Returns:\n",
    "    figure with two histograms, isi1, isi2\n",
    "\n",
    "  \"\"\"\n",
    "  plt.figure(figsize=(11, 4))\n",
    "  my_bins = np.linspace(10, 30, 20)\n",
    "  plt.subplot(121)\n",
    "  plt.hist(isi1, bins=my_bins, color='b', alpha=0.5)\n",
    "  plt.xlabel('ISI (ms)')\n",
    "  plt.ylabel('count')\n",
    "  plt.title(r'$\\sigma_{GWN}=$%.1f, CV$_{\\mathrm{isi}}$=%.3f' % (sigma1, cv1))\n",
    "\n",
    "  plt.subplot(122)\n",
    "  plt.hist(isi2, bins=my_bins, color='b', alpha=0.5)\n",
    "  plt.xlabel('ISI (ms)')\n",
    "  plt.ylabel('count')\n",
    "  plt.title(r'$\\sigma_{GWN}=$%.1f, CV$_{\\mathrm{isi}}$=%.3f' % (sigma2, cv2))\n",
    "  plt.tight_layout()\n",
    "  plt.show()\n",
    "\n",
    "\n",
    "# the function plot the raster of the Poisson spike train\n",
    "def my_raster_Poisson(range_t, spike_train, n):\n",
    "  \"\"\"\n",
    "  Generates poisson trains\n",
    "\n",
    "  Args:\n",
    "    range_t     : time sequence\n",
    "    spike_train : binary spike trains, with shape (N, Lt)\n",
    "    n           : number of Poisson trains plot\n",
    "\n",
    "  Returns:\n",
    "    Raster plot of the spike train\n",
    "  \"\"\"\n",
    "\n",
    "  # find the number of all the spike trains\n",
    "  N = spike_train.shape[0]\n",
    "\n",
    "  # n should smaller than N:\n",
    "  if n > N:\n",
    "    print('The number n exceeds the size of spike trains')\n",
    "    print('The number n is set to be the size of spike trains')\n",
    "    n = N\n",
    "\n",
    "  # plot rater\n",
    "  i = 0\n",
    "  while i < n:\n",
    "    if spike_train[i, :].sum() > 0.:\n",
    "      t_sp = range_t[spike_train[i, :] > 0.5]  # spike times\n",
    "      plt.plot(t_sp, i * np.ones(len(t_sp)), 'k|', ms=10, markeredgewidth=2)\n",
    "    i += 1\n",
    "  plt.xlim([range_t[0], range_t[-1]])\n",
    "  plt.ylim([-0.5, n + 0.5])\n",
    "  plt.xlabel('Time (ms)', fontsize=12)\n",
    "  plt.ylabel('Neuron ID', fontsize=12)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Section 1: The Leaky Integrate-and-Fire (LIF) model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:42.306202Z",
     "iopub.status.busy": "2021-05-25T01:12:42.305318Z",
     "iopub.status.idle": "2021-05-25T01:12:42.347139Z",
     "shell.execute_reply": "2021-05-25T01:12:42.347734Z"
    }
   },
   "outputs": [],
   "source": [
    "#@title Video 1: Reduced Neuron Models\n",
    "from IPython.display import YouTubeVideo\n",
    "video = YouTubeVideo(id='rSExvwCVRYg', width=854, height=480, fs=1)\n",
    "print(\"Video available at https://youtube.com/watch?v=\" + video.id)\n",
    "video"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Implementation of an LIF neuron model\n",
    "\n",
    "Now, it's your turn to implement one of the simplest mathematical model of a neuron: the leaky integrate-and-fire (LIF) model. The basic idea of LIF neuron was proposed in 1907 by Louis Édouard Lapicque, long before we understood the electrophysiology of a neuron (see a translation of [Lapicque's paper](https://pubmed.ncbi.nlm.nih.gov/17968583/) ). More details of the model can be found in the book [**Theoretical neuroscience**](http://www.gatsby.ucl.ac.uk/~dayan/book/) by Peter Dayan and Laurence F. Abbott.\n",
    "\n",
    "The subthreshold membrane potential dynamics of a LIF neuron is described by\n",
    "\n",
    "\\begin{eqnarray}\n",
    "C_m\\frac{dV}{dt} = -g_L(V-E_L) + I,\\quad (1)\n",
    "\\end{eqnarray}\n",
    "\n",
    "where $C_m$ is the membrane capacitance, $V$ is the membrane potential, $g_L$ is the leak conductance ($g_L = 1/R$, the inverse of the leak resistance $R$ mentioned in previous tutorials), $E_L$ is the resting potential, and $I$ is the external input current. \n",
    "\n",
    "Dividing both sides of the above equation by $g_L$ gives\n",
    "\n",
    "\\begin{align}\n",
    "\\tau_m\\frac{dV}{dt} = -(V-E_L) + \\frac{I}{g_L}\\,,\\quad (2)\n",
    "\\end{align}\n",
    "\n",
    "where the $\\tau_m$ is membrane time constant and is defined as $\\tau_m=C_m/g_L$. \n",
    "\n",
    "You might wonder why dividing capacitance by conductance gives units of time! Find out yourself why this is true.\n",
    "\n",
    "Below, we will use Eqn.(2) to simulate LIF neuron dynamics. \n",
    "\n",
    "If $I$ is sufficiently strong such that $V$ reaches a certain threshold value $V_{\\rm th}$, $V$ is reset to a reset potential $V_{\\rm reset}< V_{\\rm th}$, and voltage is clamped to $V_{\\rm reset}$ for $\\tau_{\\rm ref}$ ms, mimicking the refractoriness of the neuron during an action potential:\n",
    "\n",
    "\\begin{eqnarray}\n",
    "\\mathrm{if}\\quad V(t_{\\text{sp}})\\geq V_{\\rm th}&:& V(t)=V_{\\rm reset} \\text{  for } t\\in(t_{\\text{sp}}, t_{\\text{sp}} + \\tau_{\\text{ref}}]\n",
    "\\end{eqnarray}\n",
    "where $t_{\\rm sp}$ is the spike time when $V(t)$ just exceeded $V_{\\rm th}$.\n",
    "\n",
    "(__Note__: in the lecture slides, $\\theta$ corresponds to the threshold voltage $V_{th}$, and $\\Delta$ corresponds to the refractory time $\\tau_{\\rm ref}$.)\n",
    "\n",
    "Thus, the LIF model captures the facts that a neuron:\n",
    "- performs spatial and temporal integration of synaptic inputs \n",
    "- generates a spike when the voltage reaches a certain threshold\n",
    "- goes refractory during the action potential\n",
    "- has a leaky membrane \n",
    "\n",
    "The LIF model assumes that the spatial and temporal integration of inputs is linear. Also, membrane potential dynamics close to the spike threshold are much slower in LIF neurons than in real neurons. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Exercise 1: Python code to simulate the LIF neuron\n",
    "\n",
    "We now write Python code to calculate Eqn. (2) and simulate the LIF neuron dynamics. We will use the Euler method, which you saw in the linear systems case last week,  to numerically integrate Eq 2.\n",
    "\n",
    "The cell below initializes a dictionary that stores parameters of the LIF neuron model and the simulation scheme. You can use `pars=default_pars(T=simulation_time, dt=time_step)` to get the parameters (you can try to print the dictionary `pars`). Note that, `simulation_time` and `time_step` have the unit `ms`. In addition, you can add the value to a new parameter by `pars['New_param'] = value`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:42.353700Z",
     "iopub.status.busy": "2021-05-25T01:12:42.353176Z",
     "iopub.status.idle": "2021-05-25T01:12:42.356364Z",
     "shell.execute_reply": "2021-05-25T01:12:42.356765Z"
    }
   },
   "outputs": [],
   "source": [
    "# @title\n",
    "\n",
    "# @markdown Execute this code to initialize the default parameters\n",
    "\n",
    "\n",
    "def default_pars(**kwargs):\n",
    "  pars = {}\n",
    "\n",
    "  # typical neuron parameters#\n",
    "  pars['V_th'] = -55.     # spike threshold [mV]\n",
    "  pars['V_reset'] = -75.  # reset potential [mV]\n",
    "  pars['tau_m'] = 10.     # membrane time constant [ms]\n",
    "  pars['g_L'] = 10.       # leak conductance [nS]\n",
    "  pars['V_init'] = -75.   # initial potential [mV]\n",
    "  pars['E_L'] = -75.      # leak reversal potential [mV]\n",
    "  pars['tref'] = 2.       # refractory time (ms)\n",
    "\n",
    "  # simulation parameters #\n",
    "  pars['T'] = 400.  # Total duration of simulation [ms]\n",
    "  pars['dt'] = .1   # Simulation time step [ms]\n",
    "\n",
    "  # external parameters if any #\n",
    "  for k in kwargs:\n",
    "    pars[k] = kwargs[k]\n",
    "\n",
    "  pars['range_t'] = np.arange(0, pars['T'], pars['dt'])  # Vector of discretized time points [ms]\n",
    "\n",
    "  return pars\n",
    "\n",
    "\n",
    "pars = default_pars()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "The cell below defines the function to simulate the LIF neuron when receiving external current inputs. You can use `v, sp = run_LIF(pars, Iinj)` to get the membrane potential (`v`) and spike train (`sp`) given the dictionary `pars` and input current `Iinj`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "both",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:42.368269Z",
     "iopub.status.busy": "2021-05-25T01:12:42.366923Z",
     "iopub.status.idle": "2021-05-25T01:12:42.368863Z",
     "shell.execute_reply": "2021-05-25T01:12:42.369332Z"
    }
   },
   "outputs": [],
   "source": [
    "def run_LIF(pars, Iinj, stop=False):\n",
    "  \"\"\"\n",
    "  Simulate the LIF dynamics with external input current\n",
    "\n",
    "  Args:\n",
    "    pars       : parameter dictionary\n",
    "    Iinj       : input current [pA]. The injected current here can be a value\n",
    "                 or an array\n",
    "    stop       : boolean. If True, use a current pulse\n",
    "\n",
    "  Returns:\n",
    "    rec_v      : membrane potential\n",
    "    rec_sp     : spike times\n",
    "  \"\"\"\n",
    "\n",
    "  # Set parameters\n",
    "  V_th, V_reset = pars['V_th'], pars['V_reset']\n",
    "  tau_m, g_L = pars['tau_m'], pars['g_L']\n",
    "  V_init, E_L = pars['V_init'], pars['E_L']\n",
    "  dt, range_t = pars['dt'], pars['range_t']\n",
    "  Lt = range_t.size\n",
    "  tref = pars['tref']\n",
    "\n",
    "  # Initialize voltage and current\n",
    "  v = np.zeros(Lt)\n",
    "  v[0] = V_init\n",
    "  Iinj = Iinj * np.ones(Lt)\n",
    "  if stop:  # set end of current to 0 if current pulse\n",
    "    Iinj[:int(len(Iinj) / 2) - 1000] = 0\n",
    "    Iinj[int(len(Iinj) / 2) + 1000:] = 0\n",
    "  tr = 0.  # the count for refractory duration\n",
    "\n",
    "  # Simulate the LIF dynamics\n",
    "  rec_spikes = []  # record spike times\n",
    "  for it in range(Lt - 1):\n",
    "    if tr > 0:  # check for refractoriness\n",
    "      v[it] = V_reset\n",
    "      tr = tr - 1\n",
    "    elif v[it] >= V_th:  # reset voltage and record spike event\n",
    "      rec_spikes.append(it)\n",
    "      v[it] = V_reset\n",
    "      tr = tref / dt\n",
    "\n",
    "    ########################################################################\n",
    "    ## TODO for students: compute the membrane potential v, spike train sp #\n",
    "    # Fill out function and remove\n",
    "    raise NotImplementedError('Student Exercise: calculate the dv/dt and the update step!')\n",
    "    ########################################################################\n",
    "    # calculate the increment of the membrane potential\n",
    "    dv = ...\n",
    "\n",
    "    # update the membrane potential\n",
    "    v[it + 1] = ...\n",
    "\n",
    "  rec_spikes = np.array(rec_spikes) * dt\n",
    "\n",
    "  return v, rec_spikes\n",
    "\n",
    "\n",
    "pars = default_pars(T=500)\n",
    "# Uncomment below to test your function\n",
    "# v, sp = run_LIF(pars, Iinj=100, stop=True)\n",
    "# plot_volt_trace(pars, v, sp)\n",
    "# plt.ylim([-80, -60])\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:42.378515Z",
     "iopub.status.busy": "2021-05-25T01:12:42.377986Z",
     "iopub.status.idle": "2021-05-25T01:12:42.846631Z",
     "shell.execute_reply": "2021-05-25T01:12:42.846088Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove solution\n",
    "def run_LIF(pars, Iinj, stop=False):\n",
    "  \"\"\"\n",
    "  Simulate the LIF dynamics with external input current\n",
    "\n",
    "  Args:\n",
    "    pars       : parameter dictionary\n",
    "    Iinj       : input current [pA]. The injected current here can be a value\n",
    "                 or an array\n",
    "    stop       : boolean. If True, use a current pulse\n",
    "\n",
    "  Returns:\n",
    "    rec_v      : membrane potential\n",
    "    rec_sp     : spike times\n",
    "  \"\"\"\n",
    "  # Set parameters\n",
    "  V_th, V_reset = pars['V_th'], pars['V_reset']\n",
    "  tau_m, g_L = pars['tau_m'], pars['g_L']\n",
    "  V_init, E_L = pars['V_init'], pars['E_L']\n",
    "  dt, range_t = pars['dt'], pars['range_t']\n",
    "  Lt = range_t.size\n",
    "  tref = pars['tref']\n",
    "\n",
    "  # Initialize voltage and current\n",
    "  v = np.zeros(Lt)\n",
    "  v[0] = V_init\n",
    "  Iinj = Iinj * np.ones(Lt)\n",
    "  if stop:  # set end of current to 0 if current pulse\n",
    "    Iinj[:int(len(Iinj) / 2) - 1000] = 0\n",
    "    Iinj[int(len(Iinj) / 2) + 1000:] = 0\n",
    "  tr = 0.  # the count for refractory duration\n",
    "\n",
    "  # Simulate the LIF dynamics\n",
    "  rec_spikes = []  # record spike times\n",
    "  for it in range(Lt - 1):\n",
    "    if tr > 0:  # check for refractoriness\n",
    "      v[it] = V_reset\n",
    "      tr = tr - 1\n",
    "    elif v[it] >= V_th:  # reset voltage and record spike event\n",
    "      rec_spikes.append(it)\n",
    "      v[it] = V_reset\n",
    "      tr = tref / dt\n",
    "\n",
    "    # calculate the increment of the membrane potential\n",
    "    dv = (-(v[it] - E_L) + Iinj[it] / g_L) * (dt / tau_m)\n",
    "\n",
    "    # update the membrane potential\n",
    "    v[it + 1] = v[it] + dv\n",
    "\n",
    "  rec_spikes = np.array(rec_spikes) * dt\n",
    "\n",
    "  return v, rec_spikes\n",
    "\n",
    "\n",
    "pars = default_pars(T=500)\n",
    "# Uncomment below to test your function\n",
    "v, sp = run_LIF(pars, Iinj=100, stop=True)\n",
    "with plt.xkcd():\n",
    "  plot_volt_trace(pars, v, sp)\n",
    "  plt.ylim([-80, -60])\n",
    "  plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Section 2: Response of an LIF model to different types of input currents\n",
    "In the following section, we will learn how to inject direct current and white noise to study the response of an LIF neuron."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:42.852466Z",
     "iopub.status.busy": "2021-05-25T01:12:42.851943Z",
     "iopub.status.idle": "2021-05-25T01:12:42.892997Z",
     "shell.execute_reply": "2021-05-25T01:12:42.893425Z"
    }
   },
   "outputs": [],
   "source": [
    "#@title Video 2: Response of the LIF neuron to different inputs\n",
    "from IPython.display import YouTubeVideo\n",
    "video = YouTubeVideo(id='preNGdab7Kk', width=854, height=480, fs=1)\n",
    "print(\"Video available at https://youtube.com/watch?v=\" + video.id)\n",
    "video"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Section 2.1: Direct current (DC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "### LIF neuron driven by constant current\n",
    "Investigate the voltage response of the LIF neuron when receiving a DC input of 300 pA by `run_LIF` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:42.922751Z",
     "iopub.status.busy": "2021-05-25T01:12:42.914002Z",
     "iopub.status.idle": "2021-05-25T01:12:43.135763Z",
     "shell.execute_reply": "2021-05-25T01:12:43.134893Z"
    }
   },
   "outputs": [],
   "source": [
    "pars = default_pars(T=100)  # get the parameters\n",
    "\n",
    "# Run the model to obtain v and sp\n",
    "v, sp = run_LIF(pars, Iinj=300)\n",
    "\n",
    "plot_volt_trace(pars, v, sp)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "In the plot above, you see the membrane potential of an LIF neuron. You may notice that the neuron generates a spike. But this is just a cosmetic spike only for illustration purposes. In an LIF neuron, we only need to keep track of times when the neuron hit the threshold so the postsynaptic neurons can be informed of the spike. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "### Interactive Demo: Parameter exploration of DC input amplitude \n",
    "Here's an interactive demo that shows how the LIF neuron behavior changes for DC input with different amplitudes.\n",
    "\n",
    "How much DC is needed to reach the threshold (rheobase current)? How does the membrane time constant affect the frequency of the neuron?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:43.160343Z",
     "iopub.status.busy": "2021-05-25T01:12:43.159846Z",
     "iopub.status.idle": "2021-05-25T01:12:43.382413Z",
     "shell.execute_reply": "2021-05-25T01:12:43.382816Z"
    }
   },
   "outputs": [],
   "source": [
    "# @title\n",
    "# @markdown Make sure you execute this cell to enable the widget!\n",
    "my_layout.width = '450px'\n",
    "@widgets.interact(\n",
    "    I_dc=widgets.FloatSlider(50., min=0., max=300., step=10.,\n",
    "                               layout=my_layout),\n",
    "    tau_m=widgets.FloatSlider(10., min=2., max=20., step=2.,\n",
    "                                layout=my_layout)\n",
    ")\n",
    "\n",
    "def diff_DC(I_dc=200., tau_m=10.):\n",
    "  pars = default_pars(T=100.)\n",
    "  pars['tau_m'] = tau_m\n",
    "  v, sp = run_LIF(pars, Iinj=I_dc)\n",
    "  plot_volt_trace(pars, v, sp)\n",
    "  plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:43.387841Z",
     "iopub.status.busy": "2021-05-25T01:12:43.387246Z",
     "iopub.status.idle": "2021-05-25T01:12:43.390582Z",
     "shell.execute_reply": "2021-05-25T01:12:43.391007Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove explanation\n",
    "\n",
    "\"\"\"\n",
    "1. As we increase the current, we observe that at 210 pA we cross the threshold.\n",
    "\n",
    "2. As we increase the membrane time constant (slower membrane), the firing rate\n",
    "is decreased because the membrane needs more time to reach the threshold after\n",
    "the reset.\n",
    "\"\"\";"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Section 2.2: Gaussian white noise (GWN) current \n",
    "\n",
    "Given the noisy nature of neuronal activity _in vivo_, neurons usually receive complex, time-varying inputs.\n",
    "\n",
    "To mimic this, we will now investigate the neuronal response when the LIF neuron receives Gaussian white noise $\\xi(t)$ with mean \n",
    "\n",
    "\\begin{eqnarray}\n",
    "E[\\xi(t)]=\\mu=0,\n",
    "\\end{eqnarray}\n",
    "\n",
    "and autocovariance \n",
    "\\begin{eqnarray}\n",
    "E[\\xi(t)\\xi(t+\\tau)]=\\sigma_\\xi^2 \\delta(\\tau)\n",
    "\\end{eqnarray}\n",
    "\n",
    "Note that the GWN has zero mean, that is, it describes only the fluctuations of the input received by a neuron. We can thus modify our definition of GWN to have a nonzero mean value $\\mu$ that equals the DC input, since this is the average input into the cell. The cell below defines the modified gaussian white noise currents with nonzero mean $\\mu$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:43.397089Z",
     "iopub.status.busy": "2021-05-25T01:12:43.395960Z",
     "iopub.status.idle": "2021-05-25T01:12:43.397707Z",
     "shell.execute_reply": "2021-05-25T01:12:43.398154Z"
    }
   },
   "outputs": [],
   "source": [
    "#@title\n",
    "#@markdown Execute this cell to get function to generate GWN\n",
    "\n",
    "def my_GWN(pars, mu, sig, myseed=False):\n",
    "  \"\"\"\n",
    "  Function that generates Gaussian white noise input\n",
    "\n",
    "  Args:\n",
    "    pars       : parameter dictionary\n",
    "    mu         : noise baseline (mean)\n",
    "    sig        : noise amplitute (standard deviation)\n",
    "    myseed     : random seed. int or boolean\n",
    "                 the same seed will give the same\n",
    "                 random number sequence\n",
    "\n",
    "  Returns:\n",
    "    I          : Gaussian white noise input\n",
    "  \"\"\"\n",
    "\n",
    "  # Retrieve simulation parameters\n",
    "  dt, range_t = pars['dt'], pars['range_t']\n",
    "  Lt = range_t.size\n",
    "\n",
    "  # set random seed\n",
    "  # you can fix the seed of the random number generator so that the results\n",
    "  # are reliable however, when you want to generate multiple realization\n",
    "  # make sure that you change the seed for each new realization.\n",
    "  if myseed:\n",
    "      np.random.seed(seed=myseed)\n",
    "  else:\n",
    "      np.random.seed()\n",
    "\n",
    "  # generate GWN\n",
    "  # we divide here by 1000 to convert units to sec.\n",
    "  I_gwn = mu + sig * np.random.randn(Lt) / np.sqrt(dt / 1000.)\n",
    "\n",
    "  return I_gwn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "### Exercise 2: LIF neuron driven by GWN\n",
    "\n",
    "You can generate a noisy input with `my_GWN(pars, mu, sig, myseed=False)`. Here, $\\mu=250$ and $\\sigma=5$. Note that fixing the value of the random seed (e.g., `myseed=2020`) will allow you to obtain the same result every time you run this. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:43.402783Z",
     "iopub.status.busy": "2021-05-25T01:12:43.401554Z",
     "iopub.status.idle": "2021-05-25T01:12:43.403419Z",
     "shell.execute_reply": "2021-05-25T01:12:43.403922Z"
    }
   },
   "outputs": [],
   "source": [
    "pars = default_pars(T=100.)\n",
    "sig_gwn = 5.\n",
    "mu_gwn = 250.\n",
    "#######################################################\n",
    "## TODO for students: compute the noisy input I_GWN   #\n",
    "# and then compute v, sp with inject I_GWN into LIF   #\n",
    "#######################################################\n",
    "\n",
    "# Calculate the GWN current\n",
    "# I_GWN = my_GWN(pars, ..., ..., myseed=2020)\n",
    "\n",
    "# Run the model and calculate the v and the sp\n",
    "# v, sp = run_LIF(pars, ...)\n",
    "\n",
    "# Uncomment these lines after obtaining v, sp\n",
    "# plot_GWN(pars, I_GWN)\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:43.429395Z",
     "iopub.status.busy": "2021-05-25T01:12:43.407873Z",
     "iopub.status.idle": "2021-05-25T01:12:43.955186Z",
     "shell.execute_reply": "2021-05-25T01:12:43.955647Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove solution\n",
    "pars = default_pars(T=100.)\n",
    "sig_gwn = 5.\n",
    "mu_gwn = 250.\n",
    "\n",
    "# Calculate the GWN current\n",
    "I_GWN = my_GWN(pars, mu=mu_gwn, sig=sig_gwn, myseed=2020)\n",
    "\n",
    "# Run the model and calculate the v and the sp\n",
    "v, sp = run_LIF(pars, Iinj=I_GWN)\n",
    "\n",
    "with plt.xkcd():\n",
    "  plot_GWN(pars, I_GWN)\n",
    "  plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "### Interactive Demo: LIF neuron Explorer for noisy input\n",
    "\n",
    "The mean of the GWN is the amplitude of DC. Indeed, when $\\sigma = 0$, GWN is just a DC.\n",
    "\n",
    "So the question arises how does $\\sigma$ of the GWN affect the spiking behavior of the neuron. For instance we may want to know\n",
    "- how does the minimum input (i.e. $\\mu$) needed to make a neuron spike change with increase in $\\sigma$\n",
    "- how does the spike regularity change with increase in $\\sigma$\n",
    "\n",
    "To get an intuition about these questions you can use the following interactive demo that shows how the LIF neuron behavior changes for noisy input with different amplitudes (the mean $\\mu$) and fluctuation sizes ($\\sigma$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:43.981700Z",
     "iopub.status.busy": "2021-05-25T01:12:43.981118Z",
     "iopub.status.idle": "2021-05-25T01:12:44.404688Z",
     "shell.execute_reply": "2021-05-25T01:12:44.404010Z"
    }
   },
   "outputs": [],
   "source": [
    "# @title\n",
    "\n",
    "# @markdown Make sure you execute this cell to enable the widget!\n",
    "\n",
    "my_layout.width = '450px'\n",
    "@widgets.interact(\n",
    "    mu_gwn=widgets.FloatSlider(200., min=100., max=300., step=5.,\n",
    "                               layout=my_layout),\n",
    "    sig_gwn=widgets.FloatSlider(2.5, min=0., max=5., step=.5,\n",
    "                                layout=my_layout)\n",
    ")\n",
    "\n",
    "\n",
    "def diff_GWN_to_LIF(mu_gwn, sig_gwn):\n",
    "  pars = default_pars(T=100.)\n",
    "  I_GWN = my_GWN(pars, mu=mu_gwn, sig=sig_gwn)\n",
    "  v, sp = run_LIF(pars, Iinj=I_GWN)\n",
    "  plt.figure(figsize=(12, 4))\n",
    "  plt.subplot(121)\n",
    "  plt.plot(pars['range_t'][::3], I_GWN[::3], 'b')\n",
    "  plt.xlabel('Time (ms)')\n",
    "  plt.ylabel(r'$I_{GWN}$ (pA)')\n",
    "  plt.subplot(122)\n",
    "  plot_volt_trace(pars, v, sp)\n",
    "  plt.tight_layout()\n",
    "  plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:44.407461Z",
     "iopub.status.busy": "2021-05-25T01:12:44.406931Z",
     "iopub.status.idle": "2021-05-25T01:12:44.411076Z",
     "shell.execute_reply": "2021-05-25T01:12:44.411533Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove explanation\n",
    "\n",
    "\"\"\"\n",
    "If we have bigger current fluctuations (increased sigma), the minimum input needed\n",
    "to make a neuron spike is smaller as the fluctuations can help push the voltage above\n",
    "threshold.\n",
    "\n",
    "The standard deviation (or size of current fluctuations) dictates the level of\n",
    "irregularity of the spikes; the higher the sigma the more irregular the observed\n",
    "spikes.\n",
    "\n",
    "\"\"\";"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Think!\n",
    "- As we increase the input average ($\\mu$) or the input fluctuation ($\\sigma$), the spike count changes. How much can we increase the spike count, and what might be the relationship between GWN mean/std or DC value and spike count? \n",
    "\n",
    "- We have seen above that when we inject DC, the neuron spikes in a regular manner (clock like), and this regularity is reduced when GWN is injected. The question is, how irregular can we make the neurons spiking by changing the parameters of the GWN? \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Section 3: Firing rate and spike time irregularity\n",
    "When we plot the output firing rate as a function of GWN mean or DC value, it is called the input-output transfer function of the neuron (so simply F-I curve).\n",
    "\n",
    "Spike regularity can be quantified as the **coefficient of variation (CV) of the inter-spike-interval (ISI)**:\n",
    "\\begin{align}\n",
    "\\text{CV}_{\\text{ISI}} = \\frac{std(\\text{ISI})}{mean(\\text{ISI})}\n",
    "\\end{align}\n",
    "\n",
    "A Poisson train is an example of high irregularity, in which $\\textbf{CV}_{\\textbf{ISI}} \\textbf{= 1}$. And for a clocklike (regular) process we have $\\textbf{CV}_{\\textbf{ISI}} \\textbf{= 0}$ because of **std(ISI)=0**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Interactive Demo: F-I Explorer for different `sig_gwn`\n",
    "\n",
    "How does the F-I curve of the LIF neuron change as we increase the $\\sigma$ of the GWN? We can already expect that the F-I curve will be stochastic and the results will vary from one trial to another. But will there be any other change compared to the F-I curved measured using DC?\n",
    "\n",
    "Here's an interactive demo that shows how the F-I curve of a LIF neuron changes for different levels of fluctuation $\\sigma$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:44.434468Z",
     "iopub.status.busy": "2021-05-25T01:12:44.418864Z",
     "iopub.status.idle": "2021-05-25T01:12:45.468509Z",
     "shell.execute_reply": "2021-05-25T01:12:45.469018Z"
    }
   },
   "outputs": [],
   "source": [
    "# @title\n",
    "\n",
    "# @markdown Make sure you execute this cell to enable the widget!\n",
    "\n",
    "\n",
    "my_layout.width = '450px'\n",
    "@widgets.interact(\n",
    "    sig_gwn=widgets.FloatSlider(3.0, min=0., max=6., step=0.5,\n",
    "                                layout=my_layout)\n",
    ")\n",
    "\n",
    "\n",
    "def diff_std_affect_fI(sig_gwn):\n",
    "  pars = default_pars(T=1000.)\n",
    "  I_mean = np.arange(100., 400., 10.)\n",
    "  spk_count = np.zeros(len(I_mean))\n",
    "  spk_count_dc = np.zeros(len(I_mean))\n",
    "\n",
    "  for idx in range(len(I_mean)):\n",
    "      I_GWN = my_GWN(pars, mu=I_mean[idx], sig=sig_gwn, myseed=2020)\n",
    "      v, rec_spikes = run_LIF(pars, Iinj=I_GWN)\n",
    "      v_dc, rec_sp_dc = run_LIF(pars, Iinj=I_mean[idx])\n",
    "      spk_count[idx] = len(rec_spikes)\n",
    "      spk_count_dc[idx] = len(rec_sp_dc)\n",
    "\n",
    "  # Plot the F-I curve i.e. Output firing rate as a function of input mean.\n",
    "  plt.figure()\n",
    "  plt.plot(I_mean, spk_count, 'k',\n",
    "           label=r'$\\sigma_{\\mathrm{GWN}}=%.2f$' % sig_gwn)\n",
    "  plt.plot(I_mean, spk_count_dc, 'k--', alpha=0.5, lw=4, dashes=(2, 2),\n",
    "           label='DC input')\n",
    "  plt.ylabel('Spike count')\n",
    "  plt.xlabel('Average injected current (pA)')\n",
    "  plt.legend(loc='best')\n",
    "  plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:45.477177Z",
     "iopub.status.busy": "2021-05-25T01:12:45.476645Z",
     "iopub.status.idle": "2021-05-25T01:12:45.479340Z",
     "shell.execute_reply": "2021-05-25T01:12:45.479820Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove explanation\n",
    "\n",
    "\"\"\"\n",
    "Discussion: If we use a DC input, the F-I curve is deterministic, and we can\n",
    "found its shape by solving the membrane equation of the neuron. If we have GWN,\n",
    "as we increase the sigma, the F-I curve has a more linear shape, and the neuron\n",
    "reaches its threshold using less average injected current.\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "### Exercise 3: Compute $CV_{ISI}$ values\n",
    "\n",
    "As shown above, the F-I curve becomes smoother while increasing the amplitude of the fluctuation ($\\sigma$). In addition, the fluctuation can also change the irregularity of the spikes. Let's investigate the effect of $\\mu=250$ with $\\sigma=0.5$ vs $\\sigma=3$. \n",
    "\n",
    "Fill in the code below to compute ISI, then plot the histogram of the ISI and compute the $CV_{ISI}$. Note that, you can use `np.diff` to calculate ISI."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:45.486711Z",
     "iopub.status.busy": "2021-05-25T01:12:45.486133Z",
     "iopub.status.idle": "2021-05-25T01:12:45.518554Z",
     "shell.execute_reply": "2021-05-25T01:12:45.517740Z"
    }
   },
   "outputs": [],
   "source": [
    "def isi_cv_LIF(spike_times):\n",
    "  \"\"\"\n",
    "  Calculates the inter-spike intervals (isi) and\n",
    "  the coefficient of variation (cv) for a given spike_train\n",
    "\n",
    "  Args:\n",
    "    spike_times : (n, ) vector with the spike times (ndarray)\n",
    "\n",
    "  Returns:\n",
    "    isi         : (n-1,) vector with the inter-spike intervals (ms)\n",
    "    cv          : coefficient of variation of isi (float)\n",
    "\n",
    "  \"\"\"\n",
    "  ########################################################################\n",
    "  ## TODO for students: compute the membrane potential v, spike train sp #\n",
    "  # Fill out function and remove\n",
    "  raise NotImplementedError('Student Exercise: calculate the isi and the cv!')\n",
    "  ########################################################################\n",
    "  if len(spike_times) >= 2:\n",
    "    # Compute isi\n",
    "    isi = ...\n",
    "    # Compute cv\n",
    "    cv = ...\n",
    "  else:\n",
    "    isi = np.nan\n",
    "    cv = np.nan\n",
    "\n",
    "  return isi, cv\n",
    "\n",
    "\n",
    "pars = default_pars(T=1000.)\n",
    "mu_gwn = 250\n",
    "sig_gwn1 = 0.5\n",
    "sig_gwn2 = 3.0\n",
    "\n",
    "I_GWN1 = my_GWN(pars, mu=mu_gwn, sig=sig_gwn1, myseed=2020)\n",
    "_, sp1 = run_LIF(pars, Iinj=I_GWN1)\n",
    "\n",
    "I_GWN2 = my_GWN(pars, mu=mu_gwn, sig=sig_gwn2, myseed=2020)\n",
    "_, sp2 = run_LIF(pars, Iinj=I_GWN2)\n",
    "\n",
    "# Uncomment to check your function\n",
    "# isi1, cv1 = isi_cv_LIF(sp1)\n",
    "# isi2, cv2 = isi_cv_LIF(sp2)\n",
    "\n",
    "# my_hists(isi1, isi2, cv1, cv2, sig_gwn1, sig_gwn2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:45.525791Z",
     "iopub.status.busy": "2021-05-25T01:12:45.525217Z",
     "iopub.status.idle": "2021-05-25T01:12:46.230377Z",
     "shell.execute_reply": "2021-05-25T01:12:46.229928Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove solution\n",
    "def isi_cv_LIF(spike_times):\n",
    "  \"\"\"\n",
    "  Calculates the inter-spike intervals (isi) and\n",
    "  the coefficient of variation (cv) for a given spike_train\n",
    "\n",
    "  Args:\n",
    "    spike_times : (n, ) vector with the spike times (ndarray)\n",
    "\n",
    "  Returns:\n",
    "    isi         : (n-1,) vector with the inter-spike intervals (ms)\n",
    "    cv          : coefficient of variation of isi (float)\n",
    "\n",
    "  \"\"\"\n",
    "\n",
    "  if len(spike_times) >= 2:\n",
    "    # Compute isi\n",
    "    isi = np.diff(spike_times)\n",
    "    # Compute cv\n",
    "    cv = isi.std()/isi.mean()\n",
    "  else:\n",
    "    isi = np.nan\n",
    "    cv = np.nan\n",
    "\n",
    "  return isi, cv\n",
    "\n",
    "\n",
    "pars = default_pars(T=1000.)\n",
    "mu_gwn = 250\n",
    "sig_gwn1 = 0.5\n",
    "sig_gwn2 = 3.0\n",
    "\n",
    "I_GWN1 = my_GWN(pars, mu=mu_gwn, sig=sig_gwn1, myseed=2020)\n",
    "_, sp1 = run_LIF(pars, Iinj=I_GWN1)\n",
    "\n",
    "I_GWN2 = my_GWN(pars, mu=mu_gwn, sig=sig_gwn2, myseed=2020)\n",
    "_, sp2 = run_LIF(pars, Iinj=I_GWN2)\n",
    "\n",
    "# Uncomment to check your function\n",
    "isi1, cv1 = isi_cv_LIF(sp1)\n",
    "isi2, cv2 = isi_cv_LIF(sp2)\n",
    "\n",
    "with plt.xkcd():\n",
    "  my_hists(isi1, isi2, cv1, cv2, sig_gwn1, sig_gwn2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Interactive Demo: Spike irregularity explorer for different `sig_gwn`\n",
    "\n",
    "In the above illustration, we see that the CV of inter-spike-interval (ISI) distribution depends on $\\sigma$ of GWN. What about the mean of GWN, should that also affect the CV$_{\\rm ISI}$? If yes, how? Does the efficacy of $\\sigma$ in increasing the CV$_{\\rm ISI}$ depend on $\\mu$?\n",
    "\n",
    "In the following interactive demo, you will examine how different levels of fluctuation $\\sigma$ affect the CVs for different average injected currents ($\\mu$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:46.253676Z",
     "iopub.status.busy": "2021-05-25T01:12:46.252914Z",
     "iopub.status.idle": "2021-05-25T01:12:46.692231Z",
     "shell.execute_reply": "2021-05-25T01:12:46.689685Z"
    }
   },
   "outputs": [],
   "source": [
    "#@title\n",
    "\n",
    "#@markdown Make sure you execute this cell to enable the widget!\n",
    "\n",
    "my_layout.width = '450px'\n",
    "@widgets.interact(\n",
    "    sig_gwn=widgets.FloatSlider(0.0, min=0., max=10.,\n",
    "                                step=0.5, layout=my_layout)\n",
    ")\n",
    "\n",
    "\n",
    "def diff_std_affect_fI(sig_gwn):\n",
    "  pars = default_pars(T=1000.)\n",
    "  I_mean = np.arange(100., 400., 20)\n",
    "  spk_count = np.zeros(len(I_mean))\n",
    "  cv_isi = np.empty(len(I_mean))\n",
    "\n",
    "  for idx in range(len(I_mean)):\n",
    "    I_GWN = my_GWN(pars, mu=I_mean[idx], sig=sig_gwn)\n",
    "    v, rec_spikes = run_LIF(pars, Iinj=I_GWN)\n",
    "    spk_count[idx] = len(rec_spikes)\n",
    "    if len(rec_spikes) > 3:\n",
    "      isi = np.diff(rec_spikes)\n",
    "      cv_isi[idx] = np.std(isi) / np.mean(isi)\n",
    "\n",
    "  # Plot the F-I curve i.e. Output firing rate as a function of input mean.\n",
    "  plt.figure()\n",
    "  plt.plot(I_mean[spk_count > 5], cv_isi[spk_count > 5], 'bo', alpha=0.5)\n",
    "  plt.xlabel('Average injected current (pA)')\n",
    "  plt.ylabel(r'Spike irregularity ($\\mathrm{CV}_\\mathrm{ISI}$)')\n",
    "  plt.ylim(-0.1, 1.5)\n",
    "  plt.grid(True)\n",
    "  plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "### Try to answer the following:\n",
    "- Does the standard deviation of the injected current affect the F-I curve in any qualitative manner?\n",
    "- Why does increasing the mean of GWN reduce the $CV_{ISI}$?\n",
    "- If you plot spike count (or rate) vs. $CV_{ISI}$, should there be a relationship between the two? Try out yourself.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:46.697030Z",
     "iopub.status.busy": "2021-05-25T01:12:46.696447Z",
     "iopub.status.idle": "2021-05-25T01:12:46.699184Z",
     "shell.execute_reply": "2021-05-25T01:12:46.699688Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove explanation\n",
    "\n",
    "\"\"\"\n",
    "1. Yes, it does. With DC input the F-I curve has a strong non-linearity but when\n",
    "a neuron is driven with GWN, as we increase the $\\sigma$ the non-linearity is\n",
    "smoothened out. Essentially, in this case noise is acting to suppress the\n",
    "non-linearities and render a neuron as a linear system.\n",
    "\n",
    "2. (here is a short answer) When we increase the mean of the GWN, at some point\n",
    "effective input mean is above the spike threshold and then the neuron operates\n",
    "in the so called mean-driven regime -- as the input is so high all the neuron\n",
    "is does is charge up to the spike threshold and reset. This essentially gives\n",
    "almost regular spiking.\n",
    "\n",
    "3. In an LIF, high firing rates are achieved for high GWN mean. Higher the mean,\n",
    "higher the firing rate and lower the CV_ISI. So you will expect that as firing rate\n",
    "increases, spike irregularity decreases. This is because of the spike threshold.\n",
    "FOr a Poisson process there is no relationship between spike rate and spike\n",
    "irregularity.\n",
    "\"\"\";"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Section 4: Generation of Poisson type spike trains\n",
    "\n",
    "*In the next tutorials, we will often use Poisson type spike train to explore properties of neurons and synapses. Therefore, it is good to know how to generate Poisson type spike trains.*\n",
    "\n",
    "Mathematically, a spike train is a Point Process. One of the simplest models of a sequence of presynaptic pulse inputs is the Poisson process. We know that given temporal integration and refractoriness, neurons cannot behave as a Poisson Process, and Gamma Process gives a better approximation (*find out what might be the difference in the two processes*). \n",
    "Here, however, we will assume that the incoming spikes are following Poisson statistics. A question arises about how to simulate a Poisson process. The generation of the Poisson process can be realized by at least two following ways:\n",
    "\n",
    "- By definition, for a Poisson process with rate $\\lambda$, the probability of finding one event in the time window with a sufficiently small length $\\Delta t$ is $P(N = 1) = \\lambda \\Delta t$. Therefore, in each time window, we generate a uniformly distributed random variable $r \\in [0,1]$ and generate a Poisson event when $r <\\lambda \\Delta t$. This method allows us to generate Poisson distributed spikes in an online manner.\n",
    "\n",
    "\n",
    "- The interval $t_{k+1}-t_{k}$ between two Poisson events with rate $\\lambda$ follows the exponential distribution, i.e., $P(t_{k+1}-t_{k}<t) = 1 - e^{\\lambda t}$. Therefore, we only need to generate a set of exponentially distributed variables $\\{s_k\\}$ to obtain the timing of Poisson events $t_{k+1}=t_{k}+s_{k}$. In this method, we need to generate all future spikes at once. \n",
    "\n",
    "Below, we use the first method in a function `Poisson_generator`, which takes arguments `(pars, rate, n, myseed)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:46.705954Z",
     "iopub.status.busy": "2021-05-25T01:12:46.704708Z",
     "iopub.status.idle": "2021-05-25T01:12:46.706588Z",
     "shell.execute_reply": "2021-05-25T01:12:46.707041Z"
    }
   },
   "outputs": [],
   "source": [
    "# @title\n",
    "# @markdown Execute this cell to get a Poisson_generator function\n",
    "\n",
    "\n",
    "def Poisson_generator(pars, rate, n, myseed=False):\n",
    "  \"\"\"\n",
    "  Generates poisson trains\n",
    "\n",
    "  Args:\n",
    "    pars            : parameter dictionary\n",
    "    rate            : noise amplitute [Hz]\n",
    "    n               : number of Poisson trains\n",
    "    myseed          : random seed. int or boolean\n",
    "\n",
    "  Returns:\n",
    "    pre_spike_train : spike train matrix, ith row represents whether\n",
    "                      there is a spike in ith spike train over time\n",
    "                      (1 if spike, 0 otherwise)\n",
    "  \"\"\"\n",
    "\n",
    "  # Retrieve simulation parameters\n",
    "  dt, range_t = pars['dt'], pars['range_t']\n",
    "  Lt = range_t.size\n",
    "\n",
    "  # set random seed\n",
    "  if myseed:\n",
    "      np.random.seed(seed=myseed)\n",
    "  else:\n",
    "      np.random.seed()\n",
    "\n",
    "  # generate uniformly distributed random variables\n",
    "  u_rand = np.random.rand(n, Lt)\n",
    "\n",
    "  # generate Poisson train\n",
    "  poisson_train = 1. * (u_rand < rate * (dt / 1000.))\n",
    "\n",
    "  return poisson_train"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:46.711398Z",
     "iopub.status.busy": "2021-05-25T01:12:46.710846Z",
     "iopub.status.idle": "2021-05-25T01:12:47.032580Z",
     "shell.execute_reply": "2021-05-25T01:12:47.032095Z"
    }
   },
   "outputs": [],
   "source": [
    "# we can use Poisson_generator to mimic presynaptic spike trains\n",
    "pars = default_pars()\n",
    "pre_spike_train = Poisson_generator(pars, rate=10, n=100, myseed=2020)\n",
    "my_raster_Poisson(pars['range_t'], pre_spike_train, 100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "How do we make sure that the above spike trains are following Poisson statistics?\n",
    "A Poisson process must have the following properties:\n",
    "- The ratio of the mean and variance of spike count is 1\n",
    "- Inter-spike-intervals are exponentially distributed\n",
    "- Spike times are irregular i.e. $CV_{\\rm ISI} = 1$\n",
    "- Adjacent spike intervals are independent of each other.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Summary\n",
    "\n",
    "Congratulations! You've just built a leaky integrate-and-fire (LIF) neuron model from scratch, and studied its dynamics in response to various types of inputs, having:\n",
    "\n",
    "- simulated the LIF neuron model\n",
    "\n",
    "- driven the LIF neuron with external inputs, such as direct current, Gaussian white noise, and Poisson spike trains, etc.\n",
    "\n",
    "- studied how different inputs affect the  LIF neuron's output (firing rate and spike time irregularity),\n",
    "\n",
    "with a special focus on low rate and irregular firing regime to mimc real cortical neurons. The next tutorial will look at how spiking statistics may be influenced by a neuron's input statistics.\n",
    "\n",
    "However, if you have extra time, follow the section below to explore a different type of noise input.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Bonus 1: Orenstein-Uhlenbeck Process\n",
    "\n",
    "When a neuron receives spiking input, the synaptic current is Shot Noise -- which is a kind of colored noise and the spectrum of the noise determined by the synaptic kernel time constant. That is, a neuron is driven by **colored noise** and not GWN.\n",
    "\n",
    "We can model colored noise using the Ohrenstein-Uhlenbeck process - filtered white noise. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Ornstein-Uhlenbeck (OU) current\n",
    "\n",
    "We next study if the input current is temporally correlated and is modeled as an Ornstein-Uhlenbeck process $\\eta(t)$, i.e., low-pass filtered GWN with a time constant $\\tau_{\\eta}$: \n",
    "\n",
    "$$\\tau_\\eta \\frac{d}{dt}\\eta(t) = \\mu-\\eta(t) + \\sigma_\\eta\\sqrt{2\\tau_\\eta}\\xi(t).$$\n",
    "\n",
    "**Hint:** An OU process as defined above has\n",
    "\n",
    "$$E[\\eta(t)]=\\mu$$ \n",
    "\n",
    "and autocovariance \n",
    "\n",
    "$$[\\eta(t)\\eta(t+\\tau)]=\\sigma_\\eta^2e^{-|t-\\tau|/\\tau_\\eta},$$\n",
    "\n",
    "which can be used to check your code. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:47.039618Z",
     "iopub.status.busy": "2021-05-25T01:12:47.039019Z",
     "iopub.status.idle": "2021-05-25T01:12:47.042304Z",
     "shell.execute_reply": "2021-05-25T01:12:47.041877Z"
    }
   },
   "outputs": [],
   "source": [
    "# @title `my_OU(pars, mu, sig, myseed)`\n",
    "# @markdown Ececute this cell to enable the OU process\n",
    "\n",
    "\n",
    "def my_OU(pars, mu, sig, myseed=False):\n",
    "  \"\"\"\n",
    "  Function that produces Ornstein-Uhlenbeck input\n",
    "\n",
    "  Args:\n",
    "    pars       : parameter dictionary\n",
    "    sig        : noise amplitute\n",
    "    myseed     : random seed. int or boolean\n",
    "\n",
    "  Returns:\n",
    "    I_ou       : Ornstein-Uhlenbeck input current\n",
    "  \"\"\"\n",
    "\n",
    "  # Retrieve simulation parameters\n",
    "  dt, range_t = pars['dt'], pars['range_t']\n",
    "  Lt = range_t.size\n",
    "  tau_ou = pars['tau_ou']  # [ms]\n",
    "\n",
    "  # set random seed\n",
    "  if myseed:\n",
    "    np.random.seed(seed=myseed)\n",
    "  else:\n",
    "    np.random.seed()\n",
    "\n",
    "  # Initialize\n",
    "  noise = np.random.randn(Lt)\n",
    "  I_ou = np.zeros(Lt)\n",
    "  I_ou[0] = noise[0] * sig\n",
    "\n",
    "  # generate OU\n",
    "  for it in range(Lt-1):\n",
    "    I_ou[it+1] = I_ou[it] + (dt / tau_ou) * (mu - I_ou[it]) + np.sqrt(2 * dt / tau_ou) * sig * noise[it + 1]\n",
    "\n",
    "  return I_ou"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "### Interactive Demo: LIF Explorer with OU input\n",
    "In the following, we will check how a neuron responds to a noisy current that follows the statistics of an OU process."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:47.074866Z",
     "iopub.status.busy": "2021-05-25T01:12:47.074279Z",
     "iopub.status.idle": "2021-05-25T01:12:47.495907Z",
     "shell.execute_reply": "2021-05-25T01:12:47.495431Z"
    }
   },
   "outputs": [],
   "source": [
    "# @title\n",
    "\n",
    "# @markdown Remember to enable the widget by running the cell!\n",
    "\n",
    "my_layout.width = '450px'\n",
    "@widgets.interact(\n",
    "  tau_ou=widgets.FloatSlider(10.0, min=5., max=20.,\n",
    "                             step=2.5, layout=my_layout),\n",
    "  sig_ou=widgets.FloatSlider(10.0, min=5., max=40.,\n",
    "                             step=2.5, layout=my_layout),\n",
    "  mu_ou=widgets.FloatSlider(190.0, min=180., max=220.,\n",
    "                            step=2.5, layout=my_layout)\n",
    ")\n",
    "\n",
    "\n",
    "def LIF_with_OU(tau_ou=10., sig_ou=40., mu_ou=200.):\n",
    "  pars = default_pars(T=1000.)\n",
    "  pars['tau_ou'] = tau_ou  # [ms]\n",
    "\n",
    "  I_ou = my_OU(pars, mu_ou, sig_ou)\n",
    "\n",
    "  v, sp = run_LIF(pars, Iinj=I_ou)\n",
    "\n",
    "  plt.figure(figsize=(12, 4))\n",
    "  plt.subplot(121)\n",
    "  plt.plot(pars['range_t'], I_ou, 'b', lw=1.0)\n",
    "  plt.xlabel('Time (ms)')\n",
    "  plt.ylabel(r'$I_{\\mathrm{OU}}$ (pA)')\n",
    "  plt.subplot(122)\n",
    "  plot_volt_trace(pars, v, sp)\n",
    "  plt.tight_layout()\n",
    "  plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Think!\n",
    "- How does the OU type input change neuron responsiveness? \n",
    "- What do you think will happen to the spike pattern and rate if you increased or decreased the time constant of the OU process?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:47.503154Z",
     "iopub.status.busy": "2021-05-25T01:12:47.502593Z",
     "iopub.status.idle": "2021-05-25T01:12:47.507155Z",
     "shell.execute_reply": "2021-05-25T01:12:47.506563Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove explanation\n",
    "\n",
    "\"\"\"\n",
    "Discussion:\n",
    "\n",
    "In a limiting case, when the time constant of the OU process is very long and the\n",
    "input current is almost flat, we expect the firing rate to decrease and neuron\n",
    "will spike more regularly. So as the OU process time constant increases, we expect\n",
    "firing rate and CV_ISI to decrease, if all other parameters are kept constant. We\n",
    "can also relate the OU process time constant to the membrane time constant as the\n",
    "neuron membrane does the same operation. This way we can link to the very first\n",
    "interactive demo.\n",
    "\"\"\";"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Bonus 2: Generalized Integrate-and-Fire models\n",
    "\n",
    "LIF model is not the only abstraction of real neurons. If you want to learn about more realistic types of neuronal models, watch the Bonus Video!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:47.512273Z",
     "iopub.status.busy": "2021-05-25T01:12:47.511715Z",
     "iopub.status.idle": "2021-05-25T01:12:47.561544Z",
     "shell.execute_reply": "2021-05-25T01:12:47.561073Z"
    }
   },
   "outputs": [],
   "source": [
    "#@title Video 3 (Bonus): Extensions to Integrate-and-Fire models\n",
    "from IPython.display import YouTubeVideo\n",
    "video = YouTubeVideo(id='G0b6wLhuQxE', width=854, height=480, fs=1)\n",
    "print(\"Video available at https://youtube.com/watch?v=\" + video.id)\n",
    "video"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "include_colab_link": true,
   "name": "W2D3_Tutorial1",
   "provenance": [],
   "toc_visible": true
  },
  "kernel": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
